#===============================================
# (1) global parameters
#===============================================
rm(list = ls())
mort <- 0.002 # 1/d
timestep <- 1 # d
steps <- 800
#
#===============================================
# (3) life methods of the individuals
#===============================================
grow <- function(inds){
ninds <- length(inds$age)
inds$age <- inds$age + timestep
K<-rnorm(1,0.73/365,0.73/365)
Lo<-13.9-1.39+abs(rnorm(1,13.9,13.9))
tm<-k+timestep
growth<- (Lo-inds$size)*ifelse(inds$size < Lo,1, 0)*(K+K*cos(2*pi*((tm-200)/365)));#    *cos(4.28*(time-200))
inds$size <- inds$size+growth*(runif(1)-0.5)/100
inds
}
die <- function(inds) {
subset(inds,
runif(inds$age) > mort & inds$age <= 5000)
}
IBMnr<-5000
age.dist<-abs(round(rnorm(IBMnr,0,20)))
size.dist<-abs(round(rnorm(IBMnr,1,0.3),1))
#===============================================
# (4) start individuals
#===============================================
inds <- data.frame(age = age.dist,
size = size.dist)
#===============================================
# (5) life loop
#===============================================
sample.n <- NULL
sample.size <- NULL
sample.agedist <- NULL
for (k in 1:steps) {
print(paste("timestep",k))
inds <- grow(inds)
inds <- die(inds)
sample.n <- c(sample.n,list(inds$age))
sample.size <- c(sample.size,list(inds$size))
}
#===============================================
# (6) results and graphics
#===============================================
par(mfrow=c(2,1))
time <- seq(1,steps,2)
boxplot(sample.size[time],names=as.character(time), xlab="time (d)",ylab="size distribution (d)")
boxplot(sample.n[time],names=as.character(time), xlab="time (d)",ylab="age distribution (d)")
